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We demonstrate the use of a new algorithm called the 
Flat Histogram sampling algorithm for the simulation of lat- 
tice polymer systems. Thermodynamics properties, such as 
average energy or entropy and other physical quantities such 
as end-to-end distance or radius of gyration can be easily cal- 
culated using this method. Ground-state energy can also be 
determined. We also explore the accuracy and limitations of 
this method. 

02.70.Uu, 05.10.Ln, 75.40.Mg 



I. INTRODUCTION 

With the rapid rise of computing power, Monte Carlo 
methods Q have become an important tool for study- 
ing high-dimensional systems such as proteins, polymers 
and spin-glass models where many questions remain to be 
answered. While the Metropolis algorithm ||^ , due to its 
simplicity, remains the most popular choice of method, it 
faces some severe drawbacks. Firstly, the dynamics are 
slow for a class of problems that involve rugged energy 
landscapes with multiple local minima. Secondly, a se- 
ries of simulation at different temperatures is needed to 
obtain the temperature dependence of thermodynamic 
quantities. Coupled with slow dynamics, the computa- 
tion time can be prohibitive. Thirdly, it is difficult to 
calculate free energy or entropy using this method. 

While the use of histogram method ^ can alleviate 
the second problem by reweighting the canonical distri- 
bution, the accuracy is limited to a small region in a 
parameter space. Recent methods based on the direct 
computation of the density of states are capable of over- 
coming the abovementioned drawbacks of Metropolis al- 
gorithm. The multicanonical method Q is the earliest 
of these. Entropic sampling ||^ is commonly cited as 
an equivalent but simpler formulation of multicanonical 
method. The flat histogram sampling algorithm Q is 
able to generate a flat histogram in energy space simi- 
lar to multicanonical simulations, but in a simpler and 
more efhcient way. The transition matrix Monte Carlo 
method can be used to either obtain the density 

of states directly or construct the canonical distribution 
for different temperatures. The method is based upon 
the definition of a stochastic matrix, the transition ma- 
trix T{E E'). The flat histogram sampling algorithm 
is an ideal algorithm for obtaining the transition matrix 
elements through simulations. 

The transition matrix Monte Carlo method and associ- 
ated flat histogram sampling algorithm are closely related 



to the broad histogram method [|T0|,^. However, when it 
was first proposed by Oliveira et al in 1996, the dynamics 
gave incorrect results in that the method did not gener- 
ate true microcanonical averages j^]. When corrected, 
the fiip rate is identical to flat histogram's but the name 
remains a historic misnomer. A central quantity in the 
broad histogram formulation is {N{a, AE))e, the micro- 
canonical average of the number of potential moves that 
increase energy by AE = E' — E. Using this definition 
together with the requirement that moves must be re- 
versible, a general equation called the broad histogram 
equation can be derived | [T^ . 

Although the definition of {N{a, AE)) e works well for 
Ising model and can be used to construct the transition 
matrix T{E E'), this interpretation is less general 
and poses a problem for other class of problems, such as 
polymer system. We shall define a more general quantity, 
Too{E — > £"), the transition matrix at infinite temper- 
ature or simply called "infinite temperature transition 
matrix" . The matrix elements of the infinite tempera- 
ture transition matrix reduces to {N{a, AE)) e for some 
particular cases. The density of states n{E), corresponds 
to the left eigenvector of the infinite temperature transi- 
tion matrix. However, it is easier to obtain the density of 
states through another set of equations, the detailed bal- 
ance equations (analogous to the broad histogram equa- 
tion) explained later. Our procedure for solving the den- 
sity of states is also different from what is prescribed by 
the broad histogram method. 

In the following section, we shall briefly describe our 
simulation model, the HP model and its connection to 
protein folding. The subsequent section is on the tran- 
sition matrix Monte Carlo and flat histogram sampling 
algorithm. We shall discuss how it can be applied to 
polymer models. We give some numerical results in the 
fourth section and the conclusion in fifth section. 



II. THE HP MODEL 

One of the most challenging problems in computational 
biology is the problem of protein folding. Proteins are 
heteropolymer consisting of long chains of amino acids. It 
is observed that proteins generally adopt a single unique 
"native" structure. The biological function of a protein 
is often intimately dependent upon the precise geomet- 
rical structure of this folded native state. How does the 
protein encodes this unique native state in an extremely 
large conformational space? Understanding this will be 
a major breakthrough with implications in biochemistry 
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and drug design. 

It is believed that the native state hes at the global free 
energy minimum. Only a small fraction of the total con- 
formational space can be explored using high-resolution 
force-field methods. Even more computing power is re- 
quired when solvent effects are included. Hence coarse 
grained models have been developed to model the fold- 
ing process. The HP model is a commonly used 
simple lattice model with the basic assumption that it is 
the hydrophobic forces that determines the native struc- 
ture of the protein. The model recognizes only two types 
of amino acids: hydrophobic (H) and polar (P). There 
are some good arguments for supporting this assump- 
tion which we shall not elaborate Our chief con- 
cern is to study the statistical mechanical aspect of the 
model and the performance of our algorithm. Given a 
sequence of H and P, each self-avoiding chain on a two- 
dimensional lattice is counted as one conformation. Only 
non-bonded neighbouring H H contributes to the total en- 
ergy, i.e. eHH = — 1 and ehp = epp ~ 0. Under such 
conditions, low-energy conformations are compact with 
H residues residing mainly in the core and P residues on 
the outside. The principle disadvantage of this model is 
that it leads to highly degenerate ground states, espe- 
cially in three dimensions 



III. TRANSITION MATRIX MONTE CARLO 
AND FLAT HISTOGRAM SAMPLING METHOD 

In a Monte Carlo simulation, it is usual to generate 
a sequence of states {tr^, a^, . . .} using a Markov chain 
where each state denoted by a lies in the phase space of 
the model. The Markov chain is defined by a transition 
matrix W{a a'). This stochastic matrix must satisfy 
E^' W{a ^ a') = 1 and W{a a') > 0. In addition, 
we require the detailed balance condition 



P{a)W{a a') = P{cr')W{cr' -> a) 



(1) 



to guarantee an equilibrium probability distribution 
Pia), i.e. 



P{cr)W{cr-^cr')=P{a'). 



(2) 



It is useful to view the matrix W{a ^ a') as composed 
of two independent parts — selection probability S{a 
a') and acceptance rate a{<7 — *■ a'). 

W{a ^ a') = S{(J ^ (7')a{(T ^ a'). (3) 

For example, the standard Metropolis algorithm ||^ uses 

Pia') 



a{a a') = min 1, 



P{a) 



(4) 



with S{a a') usually set to a constant. We can 
easily check that it satisfies the detailed balance condi- 
tion when the selection probability S is symmetric, i.e. 



S{a a') = S{a' a). This symmetry in S can be 
relaxed for general Monte Carlo simulation [T^] but is 
needed for the flat histogram sampling algorithm. 

By summing up the detailed balance equations for all 
states a with energy E and all states a' with energy E' , 
we have 

E(a)=E E(a')=E' 

E E Wm^'-^)- (5) 

E{a')=E E{a'')=E' 

If the configuration probability distribution is a function 
of energy only, i.e. Pier) oc f{E{a)), and defining the 
transition matrix in energy as 



T{E E') = 



1 



n 



E E 



(T 



(6) 



E{cr)=E E{a')=E' 



we have 

n{E)f{E)T{E ^ E') = n{E')f{E')T{E' ^ E). (7) 

The transition matrix T{E E') is also a stochastic 
matrix with the histogram h{E) = n{E)f{E)/Z being 
the equilibrium distribution: 



E 



h{E)T{E ^ E') = h{E'). 



(8) 



Since the acceptance rate a(cr a') in Eq. (g) is the 
same for configurations with a fixed energy, T(E — > E') 
can also be written as a product of two independent 
factors — the infinite temperature transition matrix 
Toc{E E') and acceptance rate in energy a{E E'): 

T{E^E')=T^{E^E')a{E-^E'), E ^ E' , (9) 

where a{E — > E') can be any acceptance rate satisfying 
Eq. (|I|), the detailed balance equation, e.g. Metropolis 
acceptance rate min(l, /(_£')//(£')), and 



T^{E^E') = 



E E 

E{a)=E E{<y')=E' 



S{a a') 
n{E) 



(10) 



When S{(J a') is taken to be a constant, the ex- 
pression can be simplified. For example, in spin systems 
we usually pick a spin randomly to be flipped so that 
S{a — > <t') — 1/N, for a and a' related by a single spin 
flip and zero otherwise, where N is the total number of 
spins. In this case 



T^{E^E') 



i{E) 



E 



N{a,AE) {N{a,AE))E 



E{a)=E 



N 



N 



(11) 



where AE = E' — E and N{a, AE) represents the num- 
ber of spins that changes the energy of the current con- 
figuration a by AE when fiipped, i.e. N{a, AE)/N = 
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Substituting the Eq. (|) into Eq. (0), and using the 
equation 

f{E)a{E^E')^fiE')a{E'~.E), (12) 

which is derived from the detailed balance equation 
Eq. (|l]) and Eq. (||) with the requirement that S{a — > a') 
is symmetric, which also implies that moves must be re- 
versible, we obtained a set of equations 

n{E)TUE^ E')^n{E')Too{E' ^ E). (13) 

When 5((T — > (t') is fixed to a constant, we can write 

n{E){N{a, AE))e = n{E'){Ni<j', ~AE))e'. (14) 

This equation, known as broad histogram equation, was 
first presented by Oliveira in [ll| ] and by Berg and Hans- 
mann jlSj. Their derivation is based on property that 
moves between configurations are reversible and is some- 
what simpler than our arguments above. 

However, the interpretation of N{a, AE) cannot be ap- 
plied when the selection probability is not a constant as 
the definition of N{a, AE) limits it to be an integer. The 
general definition of Too{E E') that can be applied to 
any choice oi S{a ^ a') is given by Eq. (p^. Eq. ( p^ is a 
general equation with two important assumptions made 
in our formulation: the probability of every configura- 
tion is a function of its energy only and allowable moves 
between any pair of configurations are selected with the 
same probability in both directions. 

The energy detailed balance equation can thus also be 
written as 



h{E)T^{E ^ E')a{E ^ E') = 

h{E')T^{E'^ E)a{E' 



E). (15) 



If we require that the histogram is constant or flat, i.e. 
h{E) = h{E'), we choose the acceptance rate 



a{E -> E') = min 1 



T^jE' ^ E) 
T^{E-^E') 



(16) 



This equation was first proposed in p] and ||]. While 
equivalent to the entropic sampling [g[ acceptance rate 
min(l, n(£')/n(£^')), the simulation procedure is differ- 
ent. We use a cumulative average to construct the ac- 
ceptance rate during simulation, i.e. 



Too{E 



1 



E{cr''}=E' 



(17) 



where H(E) = J2T=i^E{ai),E is the cumulative his- 
togram for a given energy, m is the total number of sam- 
ples generated so far, cr* is the configuration at step i 
of the simulation and ct'* are the available configurations 
that can be reached from cr* in one move. Whenever data 
is unavailable to compute the acceptance rate, we sim- 
ply accept the move in order to sample the unexplored 
states. 




FIG. 1. Three types of moves, 
crankshaft move. 



i.e. end, corner and 



To model protein folding by a Markov process, it is 
necessary to first define the move set. There are sev- 
eral choices and it has been shown that different move 
sets can affect kinetic quantities, like the mean first pas- 
sage time [|l9|. We use a local move set consisting of 
end, corner and the crankshaft move [ pOpl| ]. These are 
shown in Fig. |l]. Note that the end moves are restricted 
to 90° rotations and thus can have one or two possible 
positions depending on current configuration. The set 
of valid moves from current conformation a to new con- 
formation a' is those that can be performed using end, 
corner or crankshaft move, while preserving the excluded 
volume constraint. 

The choice of the move set used will affect the sampling 
space of configurations and also the correlation time of 
our algorithm. It can be shown that our set of local 
moves is non-ergodic, meaning that it does not generate 
all possible self-avoiding chains on the lattice. However, 
the number of such configurations is probably negligi- 
ble compared to our sampling error pl| . We can also 
consider the definition of our model to include only those 
states accessible by the move set specified. As long as our 
native state is accessible, there will not be a problem. 

The ergodicity problem can be overcome by pivot 
moves |22[ . Pivot moves are defined by randomly choos- 
ing a monomer as a pivot and rotating or reflecting one 
segment of the chain with the pivot as origin. However, 
pivot moves also fills more entry in our transition matrix 
making it less banded. Moreover, it is also found that 
pivot moves do not necessary lead to faster dynamics in 
the simulation |15|. 

Although some studies only consider conformations as 
distinct when they are not related by reflection or rota- 
tional symmetry, the choice should not matter since they 
differ by a factor of 8 in counting the number of states. 
For the special configurations on a straight line, the fac- 
tor is 4. However, such configurations have the highest 
energy and are thus almost negligible. 

While the selection probability is commonly taken to 
be 1/A^ in spin systems, we have some freedom in decid- 
ing the specific selection probabilities based on different 
moves in polymer systems. We list three possibilities. 

• Given a configuration, we can construct a list of all 
possible valid moves satisfying the excluded volume 
constraint. Each move is selected with a fixed equal 
probability. In practice, we assign moves into a list 
that can contain up to M moves. Since M is the up- 
per bound to maximum possible valid moves avail- 
able to any configuration, the remaining unassigned 
moves are considered invalid. A move is picked at 
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random. If it is valid, we accept the move using the 
acceptance ratio in Eq. (|l6|). In our simulation, we 
set M equal to N , the number of monomers. 

• Given a configuration, we pick a monomer at ran- 
dom. If it is an end monomer, it is able to ro- 
tate to two possible positions. [] We select one ran- 
domly without considering the excluded volume 
constraint. If it is not an end monomer, then corner 
or crankshaft moves are possible. These two moves 
are mutually exclusive. Once a move is selected, we 
check the excluded volume constraint. If a move 
results in moving a monomer into an already oc- 
cupied position, the move is rejected. Otherwise, 
we accept the move using the acceptance ratio in 
Eq. ( 16 ) . For cases where no moves are possible for 
the monomer picked, the current configuration is 
included in the average and the time step is incre- 
mented by one. This method of selecting move is 
used in Ref. |^^. It relies on the fact that in mov- 
ing to a new configuration using a particular type 
of move, only the same type of move in reverse can 
bring it back to the original configuration. Thus, 
the move set must be chosen with this property to 
ensure that selection probability is symmetric. 

• We can designate end moves and corner moves as J- 
monomei move and crankshaft move as 2-monomer 
move. With probability e.g. 0.2, we choose to 
perform 1-monomer move; if an end monomer is 
picked, end moves are selected and corner moves 
otherwise, in the same manner as above (see foot- 
note). With remaining probability 0.8, pick a 
monomer from 1 to iV — 3, a 2-monomer move is se- 
lected if monomers from itoi+3 forms a crankshaft 
otherwise a move is unsuccessful. This method is 
used in Ref. ||. 

For compact configurations, the rejection rate for the 
second method would be very high. We use the first 
method, although a large M can also make the method 
inefficient for compact (thus low-energy) configurations 
since the number of possible moves is low compared to 
the value of M. However, this can be overcome if an 
N-fold way simulation is done. A move is always 
accepted in the N-fold way and the average lifetime of a 
configuration is taken into account when averaging. The 
first method also has the advantage of saving some com- 
putations. Since the selection probability is a constant, 
we can just add a constant (equivalent to counting moves) 
when calculating Too{E E') in Eq. The hst of 

moves can also be used in constructing an N-fold way 



simulation. For other choices of selection probability, we 
will have to calculate S{a a') explicitly for each move 
before adding. Once Toc{E — > £") is sampled, we solve 
Eq. (|l3|). Since n{E) varies by a huge order of magnitude, 
we solved for lnn(£') instead. Broad histogram method 
uses a forward difference scheme of integration. We solve 
it instead using a least squares method. When multiple 
simulations are performed, we can view it as an optimiza- 
tion problem taking the variance of sampling data into 
account M. 



IV. NUMERICAL RESULTS 

We present results on a sequence with 14 monomers 
with the sequence HHHPHPHPPHPHPH. Through enu- 
meration, we found that this sequence have a unique 
ground state (i.e. 8 possible configurations in our count- 
ing) of energy —7. The full density of states is given in 
Table |. 




FIG. 2. Native state for the sequence HHHPHPHPPHPHPH. 



TABLE I. Density of states for the se- 
quence HHHPHPHPPHPHPH. ii{E) is Monte Carlo results 
using flat histogram sampling algorithm and n{E) is through 
enumeration. 

% error 
7M 
5.00 
0.89 
1.56 
1.38 
4.73 
1.42 




E 


n{E) 


h{E) 





581340 


540416.37 


-1 


228416 


217016.11 


-2 


56344 


55837.55 


-3 


12472 


12666.23 


-4 


2432 


2465.45 


-5 


464 


485.93 


-6 


24 


23.66 


-7 


8 


8 



'^This is different from our case where 180° rotations are not 
allowed. It is necessary for the number of possible moves of 
each type be unambiguous for this case. 
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FIG. 3. The average energy per monomer against temper- 
ature. The insert shows the relative error between using enu- 
merated n{E) and each method. 



The average energy per monomer and radius of gyra- 
tion are plotted in Fig. ^ and ^ and are compared with 
the single histogram method and the Metropolis algo- 
rithm. We used 10^ Monte Carlo steps for the flat his- 
togram sampling algorithm and each temperature point 
of Metropolis algorithm and also the single histogram 
reweighted at T = 0.75. There are about 140 tempera- 
ture points in the Metropolis simulation, which requires 
around a 100 times more computing time compared to 
the flat histogram simulation. 
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FIG. 4. The radius of gyration against temperature. The 
insert shows the relative error between using enumerated n{E) 
and each method. 



The plots indicate that the Metropolis algorithm be- 
comes unreliable below around T = 0.5. This slowing 
down in dynamics is also found in (23| and attributed to 



the increasingly deep kinetic traps with decreasing tem- 
perature. The flat histogram sampling algorithm is unaf- 
fected by this effect with slightly better accuracy for low 
temperature range. The single histogram method also 
produces roughly the same degree of accuracy. Here we 
do not observe the reweighting errors due to the expo- 
nential decay of the canonical distribution because the 
energy spectrum is narrow and thus adequately sampled. 
The single histogram method works well only in such a 
situation. The entropy which can be easily calculated 
from the knowledge of n{E), is shown in Fig. ||. It is 
difficult to obtain this from Metropolis simulation. 
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FIG. 5. The average entropy per monomer against temper- 
ature. The insert shows the relative error. 



Finding the energy of the native state is also an impor- 
tant task in protein folding. To generate the native states 
using Metropolis algorithm, the temperature must be low 
enough so that the canonical distribution covers the low 
energy range adequately. However, the canonical distri- 
bution has width \/iV and low temperature simulation 
causes the system to be trapped in local minima. Vari- 
ous methods, such as genetic algorithm |^ and methods 
employing heuristic ||27|| have been proposed to overcome 
this problem. Often, an annealing schedule is adopted 
whereby the temperature is lowered as the simulation 
progress. There is no standard way and considerable 
trial-and-error is necessary. 

The flat histogram sampling algorithm can be used as 
a method for determining the native state energy. Since 
the energy barriers no longer exist, we expect a random 
walk along the energy scale in the ideal case. An advan- 
tage is that in most polymer models, the energy range 
do not increase rapidly with system size. We also do not 
have to devise any annealing schedule or adjust many 
parameters. We can therefore use our algorithm for de- 
termining the native state energy. Although we cannot 
attach any physical signiflcance to the time for flnding the 



5 



native state since the ensemble is non-Boltzmann (multi- 
canonical) , it is still useful for analyzing the performance 
of our algorithm. This native state time tq, the time to 
reach the native state from an unfolded conformation, is 
shown for four sequences in Table ^ We select the first 
two sequence to have unique native state. The other two 
sequences were taken from [ p6[ . 

TABLE II. Sequence of HP monomers used in simulation 
with their lowest energy and native state time averaged from 
10* simulations. The last column is the standard deviation of 

TO- 



N 


Sequence 


Eo 


To 




10 


HPHPPHPPHH 


-4 


339.7 


405.7 


14 


HHHPHPHPPHPHPH 


-7 


5641.4 


6340.4 


20 


HPHPPHHPHPPHPHHPPHPH 


-9 


81788.1 


85478.5 


25 


PPHPPHHPPPPHHPPPPHHPPPPHH 


-8 


196137.9 


318722.3 



The tunnelling time can also be used as a measure of 
the efhciency of our algorithm. We denote t„, the "up" 
tunnelling time as the average MCS taken for a state with 
minimum energy to reach a state with maximum energy 
while Td, the "down" tunnelling time is for the opposite 
direction. These are shown in Table [11. Unlike spin 



systems, where the two tunnelling times are the same 
due to the symmetry in the Hamiltonian, it is faster to 
tunnel to higher energies than to lower energies. 

Fig. H shows the general behaviour of the native state 
time and tunnelling times as the size increases. It is 
usual in disorder systems to average over different ran- 
dom realizations corresponding to different sets of cou- 
pling constants or sequences. However, it is recognized 
that proteins are not random sequences since they fold 
into unique native states with specific properties. The 
probability of selecting a sequence with this property 
from all possible sequences is very small. This leads 
to the problem of designing sequences with protein-like 
properties. 

TABLE III. Tunnelling times, Tu and Td for different se- 
quences. 



iV 


Eo 


Tu 


std 


count 


Td 


std 


count 


10 


-4 


66.8 


60.7 


4272 


167.3 


208.3 


4271 


14 


-7 


592.6 


551.5 


2615 


3229.3 


3940.6 


2615 


20 


-9 


2224.3 


3612.3 


638 


13427.8 


13976.7 


638 


25 


-8 


8418.0 


9616.3 


190 


44214.1 


61316.7 


190 



10° 
10^ 
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10= 
10* 
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10^ 
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10" 



Oz, 
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-- z~n' 
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100 



FIG. 6. The points are native state time and tunnelling 
times obtained from simulation. The straight lines are fits to 
a power law r oc A'^*'. 



We fit the times to find the minimum energy (native) 
states and the tunnelling times for the four sequences to 
a power law. The time to find native state follows ap- 
proximately tq oc 7V^'^. The "up" and "down" tunnelling 
time gives the fitted parameter t„ oc N^'^ and oc N^^^ 
respectively. This suggests that the algorithm takes in- 
creasingly longer time to reach the lowest energy states 
but moves easily towards the upper energy levels. It re- 
flects our observation that the flat histogram sampling 
algorithm does not scale very well for longer chains es- 
pecially when the density of states increase sharply with 
energy. This also implies that the performance is not 
ideal for sequences with a unique native ground state. 
We note that our simulation is non-Markovian and thus 
the convergence is difficult to analyze. This can lead 
to detailed balance violation |2^. However, this problem 
can be alleviated by a two pass simulation. The first pass 
is the same as before. The second pass uses a fixed flip 
rate, mm{l, n{E)/n{E')), obtained from the first pass. 



V. CONCLUSION 

We have shown that the flat histogram sampling al- 
gorithm, which was flrst proposed and implemented for 
spin systems, can be used for the simulation of lattice 
polymer systems. We give some measures of its accuracy 
and also efficiency in terms of thermodynamics proper- 
ties, native state time and tunnelling times. The current 
implementation is useful for up to about 20 monomer 
HP chain. However, we would like to emphasize that the 
flat histogram sampling algorithm is still vastly superior 
to Metropolis algorithm, especially in terms of accuracy 
for low temperature properties. The simulation time is 
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modest as there is no need for simulation at each tem- 
perature point. It also has the advantage of easily ob- 
taining the density of states for free energy or entropy 
calculations. While the algorithm is rather basic, there 
are improvements to be made such as the extensions and 
modifications proposed in Ref. 
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